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Ab str ac t 

This article describes the application of the Spectrum™ 
Solver in Multidisciplinary Analysis (MDA). Spectrum, 
a multiphysics simulation software based on the finite 
element method, addresses compressible and incom- 
pressible fluid flow, structural, and thermal modeling as 
well as the interaction between these disciplines. Multi- 
physics simulation is based on a single computational 
framework for the modeling of multiple interacting 
physical phenomena. Interaction constraints are 
enforced in a fully-coupled manner using the aug- 
mented-Lagrangian method. Within the multiphysics 
framework, the finite element treatment of fluids is 
based on the Galerkin-Least-Squares (GLS) method 
with discontinuity capturing operators. The arbitrary- 
Lagrangian-Eulerian method is utilized to account for 
deformable fluid domains. The finite element treatment 
of solids and structures is based on the Hu-Washizu 
variational principle. The multiphysics architecture 
lends itself naturally to high-performance parallel com- 
puting. Aeroelastic, propulsion, thermal management 
and manufacturing applications are presented. 

1 . Introduction 

Many of the current directions in product design and 
analysis are driven by competitive and regulatory con- 
straints, such as the need to shorten design cycles, 
reduce cost, meet increasingly stringent government 
regulations, improve quality and safety, and reduce 
environmental impact. These directions have increased 
the need for accurate product and component simulation 
and pressed analysts for simulations of unprecedented 
scale and complexity. For example, aeroelastic model- 
ing (figure 1), involving strong coupling of fluids and 
structures, is an essential element in the design process 
of aircraft. Thermal management issues arising from 
automotive, aeropropulsion and aerospace design and 
manufacturing applications require the coupled analysis 
of fluids, thermal and structural models. In general, 
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flow-induced vibrations are a primary concern in the 
aerospace, automotive and defense industries. 



Figure 1* Aeroelastic simulation of the ARW2 
wing configuration. 

The high-fidelity simulation of multidisciplinary inter- 
actions is computationally expensive. While single pro- 
cessor computers improve every year, the quantum 
increase in computational needs requires a quantum 
increase in computational power. The hardware industry 
has answered this increased need with affordable high- 
performance parallel architectures which are built on 
commodity parts and compatible with workstation sys- 
tems. The multiphysics architecture used in this work 
lends itself naturally to high-performance parallel com- 
puting. Coarse grain parallel processing is utilized 
through the SPMD paradigm with domain decomposi- 
tion. Within each problem subdomain, super-scalar, vec- 
tor processing and cache efficiency is leveraged with 
element blocking schemes. 

2. The Multiphvsics Problem Model 

The multiphysics architecture supports multiple physi- 
cal interactions through a data model defined as a hierar- 
chical tree of regions and interfaces 1 . Regions of a 
problem are used to separate the different physics being 
analyzed over the spatial domain. For example, in a 
fluid-structure interaction (FSI) problem, the fluid 
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domain is one region and the solid structure is another 
(figure 2). Interfaces are used to enforce the coupling 
constraints between the different regions. 

Fluid Region Interface Solid Region 


Figure 2. Multiphysics problem domain. 


ALE boundary conditions ensure that the deforming 
mesh conforms to both the stationary and moving 
boundaries. Within the interior of the fluid domain, 
mesh movement is modeled by the equations of large 
deformation elasticity. In effect, this model computes 
the position of the interior nodes as if all nearest nodal 
neighbors were coupled by an elastic medium and sets 
the positions and velocities of the boundary nodes to 
exactly match the boundary motions. Note that this 
model is a purely mathematical construct; it is a method 
for updating the mesh in a way that has a reasonable 
chance of maintaining mesh integrity. Thus, terms such 
as Cauchy stresses, elastic moduli, and so forth do not 
have the usual physical interpretation within this con- 
text. 




Because of the varying discretization requirements of 
the different physical phenomena, this approach is 
designed to allow variable mesh densities and element 
topologies at the region interfaces. The finite element 
formulation admits (automatically generated) unstruc- 
tured meshes. For example, in a typical FSI problem, a 
highly refined tetrahedral mesh may be used to model 
the fluid domain with a relatively coarse hexahedral dis- 
cretization representing the solid region. These discreti- 
zations do not, in general, coincide at the shared region 
boundaries. 


To take mesh movement into account, the convective or 
Euler flux in equation (1) is replaced by the ALE con- 
vective flux. The ALE convective flux is related to the 
Euler flux by the equation 

fALE = p *onv_ u ALe u (2) 

ALE 

where u is the velocity field of the mesh. The Euler 
Jacobians of equation (1) become 


- ALE 

A; 


A ALE A 

= A i ~ U: Aq 


(3) 


3 a Fluid Re gi ons 

The finite element treatment of fluid regions within this 
framework is based on the Galerkin-Least-Squares 
(GLS) method with discontinuity capturing operators 2 3 " 5 . 
Reynolds-averaged and Large-Eddy-Simulation models 
are utilized for turbulence simulation 6 . In this treatment, 
the compressible flow formulation makes use of physi- 
cal entropy variables, V . With these variables, the fluid 
conservation laws are expressed in symmetric form 
which intrinsically expresses the mathematical and 
physical stability provided by the second law of thermo- 
dynamics. Specifically, the Navier-Stokes equation sys- 
tem is expressed in the form 

— — — c rr 

VW.< * <» 

In turn, the finite element techniques employed herein 
inherit this fundamental stability and convergence 
proofs are available 7 . 

Multiphysics problems often require the movement of 
the computational fluid domain in response to the defor- 
mation of the common solid region boundaries. The 
arbitrary-Lagrangian-Eulerian (ALE) method is utilized 
to account for the deformations in fluid domains 8 . 


4. Solid Regions 

The finite element treatment of solid regions within this 
framework employs a 3-field formulation based on the 
Hu-Washizu variational principle. This method is well 
documented in the literature to address numerical lock- 
ing phenomena 10 . The kinematic description admits 
small and finite deformations and strains. Linear and 
nonlinear material models are used for the constitutive 
relations with thermo-mechanical coupling. 

The balance of linear momentum in a solid continuum 
can be expressed for the current and reference configu- 
rations in terms of the Cauchy stress, a and the first 
Piola-Kirchhoff stress, P , respectively. 

diva ■+■ pb = pv 

m (4) 

DivP * Po b m = Po* 

where p is the mass density, b m is the body force, v is 
the part cle velocity and the 0 subscript denotes a quan- 
tity in the reference configuration. 

In orde r to address incompressibility locking, a mixed 
method which modifies the interpolation of the defor- 
mation gradient is used 10 . The modified deformation 
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gradient is based on a separation of the deformation gra- 
dient F into volumetric and deviatoric parts. 


constraints are enforced with the augmented Lagrangian 
formulation. 


F 


= F vol F dev 


( 5 ) 


In constructing the modified deformation gradient, a 
mixed treatment replaces the volumetric part. The modi- 
fied deformation tensor, F , is defined using a new field 
variable, 0 , for the determinant of the deformation gra- 
dient. 

The structural formulations (beams and shells) are 
expressed in resultant form 3 3 ' 32 . The exponential map is 
employed for rotational updates which are geometri- 
cally exact and singularity free. The structural elements 
are coupled to general material models using numerical 
integration of the constitutive relations through the 
thickness direction. 


5. Heat Transfer 

The energy balance equation for heat transfer in a solid 
body has the form 


dT 

P c v -fo = -d'v? + r 


(6) 



Figure 3. Mesh discretization for the aeroelastic 
simulation of the ARW2 wing. The fluid and solid 
surface grids vary in topology and density. 

7. Parallel Processing 


For isotropic heat conduction using Fourier’s Law in a 
solid material with constant properties, this reduces to 

= k div(gradT) +r (7) 

Density, p , specific heat, c v , and thermal conductivity, 
k , are material-dependent variables. The origin of the 
volumetric heat source, r , differs depending on the type 
of analysis, e.g., in a thermomechanical problem, the 
heat source can take the form of either heat dissipation 
from inelastic deformations of the body, or structural 
heating from thermal strains and temperature variations 
in the material properties. 

Convective heat flux boundary conditions are applied to 
element surfaces on a boundary T ^ . A convective heat 
flux describes the heat flow from a solid body to a sur- 
rounding fluid using Newton’s Law of cooling. 

6. Multiphvsics Interfaces 

Interactions between regions are enforced with unstruc- 
tured mesh interfaces (figure 2). This approach is 
designed to allow variable mesh densities and element 
topologies at the region interfaces (figure 3). 

A slave-master algorithm is used to define the discrete 
interface constraints of multiphysics problems 9 . These 


The same architecture that supports multiple physics 
simulation naturally and cleanly supports independent 
and parallel computation. For each problem, we main- 
tain the concept of multiple subdomains. A subdomain 
is a collection of regions which are uniform with respect 
to linear solution technology (e.g., iterative, direct, etc.). 



Figure 4. Scalability of a 250K element socket 
flow. Near- theoretical scaling with superlinearity at 
12 and 16 subdomains is observed. 


The Same-Program Multiple Data (SPMD) model is 
implemented via message-processing. In this case, there 
is support for a wide range of hardware platforms with 
the standard message-passing interfaces of PVM and 
MPI. This programming model has the advantage of 
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allowing the reuse of large portions of code from the 
uniprocessor version. Excellent performance has been 
demonstrated on a number of different parallel systems 
for large, multiphysics simulations (e.g., figure 4). The 
relative ease of programming and support as well as the 
portability, scalability and wide availability of systems 
supporting SPMD make it the preferred model for our 
applications. 



Figure 5. Domain decomposition illustration. 

The decomposition of the computational domain (figure 
5) is an active area of research which has produced 
many algorithms and general purpose tools 13 " 15 . Many 
decomposition methods commonly use the recursive 
spectral bisection (RSB) approach 13 . A significant cost 
of the RSB method is associated with the computation 
of eigenvectors of a Laplacian matrix constructed from 
the adjacency structure of the mesh. Hendrickson and 
Leland introduced a multi-level implementation for the 
construction of the Laplacian matrix, resulting in signif- 
icant CPU performance improvement 14 . Karypis and 
Kumar present a rigorous analysis of multilevel meth- 
ods and demonstrate analytically their effectiveness 15 . 

A common feature of most domain decomposition 
research is the focus on single homogeneous grid appli- 
cations. In order to support multiphysics simulations, 
algorithms which extend the multilevel method of Kary- 
pis and Kumar to heterogeneous interfaced discretiza- 
tions are used in the present approach. 

To take full advantage of the architecture described 
above, a multi-subdomain solver is employed to solve 
the matrix set of equations which result from the dis- 
crete finite element problem. In this context, the solver 
is composed of one global solver and a set of local sub- 
domain solvers. At the subdomain level, the local solver 
may be explicit, implicit iterative or implicit direct. The 
local subdomain solver may vary from subdomain to 
subdomain. Hie global solver must be implicit iterative. 

Two iterative solvers are used for multiphysics prob- 
lems: the preconditioned conjugate gradient (CG) and 
generalized minimum residual (GMRES) methods. CG 
is used for symmetric systems which arise, e.g., from 


solid linear momentum, heat transfer and mesh move- 
ment equations. GMRES is used for non-symmetric sys- 
tems which arise, e.g., from the fluid linear momentum, 
thermal, scalar transport, and turbulence equations. 


& Multiphvsics Applications 

Four applications are presented here to demonstrate the 
applicability of multiphysics simulation in thermal man- 
agement and fluid-structure interaction problems. 

8. 1 . Panel Flutter 


A 2-D panel flutter problem is used to validate the meth- 
ods developed in this work and illustrate the difference 
between linear and nonlinear aeroelastic computations. 
The problem consists of supersonic flow over a flat plate 
which is clamped on both ends (figure 6). The panel is 
initially flat, with a constant freestream pressure on both 
sides. An analytical solution exists for the linear 
aeroelasticity problem: with the panel characteristics 
listed in the table below, the stability limit (i.e, the flow 
speed above which flutter occurs) for the panel is at 
Mach 2. 


Property 

Value 

Property 

Value 

Young’s Modulus 

77.28 GPa 

Length 

0.5 m 

Poisson’s Ratio 

0.33 

Thickness 

1 .35 mm 

Density 

2,710 kg/m 3 

Pressure 

28 kPa 

^ A A A/N/v A /\AAAA a a 1 






Figure 6. Mesh discretization of the panel flutter 
problem. 

To perurb the instability, the lower side pressure is 
dropped by 0.1% for 4 milliseconds. Thereafter, the 
pressure is brought back to 28 kPa. The flow is assumed 
inviscid (compressible Euler equations). Within each 
time-step, the FSI equations are solved in a staggered 
order. The multi-subdomain solver is used for this anal- 
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ysis. The fluid region is one subdomain with the solid 
and interface regions combined in a second subdomain. 
In the fluid subdomain, GMRES is used for the flow 
equations, and conjugate gradients are used for the mesh 
movement equations. A direct solver is used on the sec- 
ond subdomain. 



Figure 7. Time histories of the plate deformation at 
different flow speeds. 


Figure 7 plots the time history of the displacement of the 
panel for several flow speeds using linear elasticity. As 
seen in the plots, the oscillations of the plate for Mach 
numbers 1.90 and 1.95 are damped, while the oscilla- 
tions for Mach numbers 2.00 and 2.05 grow unbounded. 
The case of Mach 1 .98 appears critically damped. 



Figure 8. Comparison of the linear and nonlinear 
aeroelastic response of the panel at M=2.3. Linear 
theory predicts non-physical exponential growth. 

The difference between the linear and nonlinear 
aeroelastic solutions is depicted in figure 8 for a Mach 
number of 2.3. The linear solution (consistent with the 
analytical solution) predicts an unbounded exponential 
growth of the oscillations of the panel at this unstable 
speed. The nonlinear solution, however, shows a limit 
cycle (a result often observed in experiments). The dif- 
ference between the two solutions can be attributed to 
the coupling between membrane and bending stresses in 


the panel with nonlinear elasticity. This phenomenon 
increases the effective stiffness of the panel with its 
deformation, which changes the behavior of the system. 

Many aerospace components (particularly in military 
applications) are used beyond the stability limits pre- 
dicted by the linear theory. The methods developed in 
this work offer a tool which takes advantage of high- 
performance computing to predict limit-cycle behavior 
with nonlinear aeroe las ti city. 

8. 2, A GARD 44 $.$ Win g Flutter 

The AGARD 445.6 standard aeroelastic configuration 16 
is used to assess the 3-dimensional application of this 
method for flutter prediction. The wing is a half-span 
wind tunnel wall mounted model with a quarter-chord 
sweep angle of 45 degrees. The configuration has a 
panel aspect ratio of l .65 and a tapper ratio of 0.66. The 
test model is made of mahogany with holes drilled in 
and filled with foam to add flexibility. 



Figure 9 . First mode of vibration (9.6 Hz) for the 
AGARD 445.6 wing. 



Figure 10. Close-up of the cross-section at the 
symmetry plane of the fluid-structure mesh. 

The first mode of vibration of the wing is depicted in 
figure 9. A mesh consisting of 1 18,480 tetrahedral ele- 
ments and 22,014 nodes is used to discretize the flow 
domain around the wing. The panel is modeled with 520 
hexahedral elements and 770 nodes. A cross-section of 
the mesh at the symmetry plane is shown in figure 10. 
The wing is perturbed with a combination of its first two 
vibration modes applied as initial velocities. No struc- 
tural damping is used. The oscillations of the structure 
are then used to compute the damping coefficients. Fig- 
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ure 1 1 shows the deformation of the wing (magnified) at 
two different instances during the oscillation. Also 
shown are contour lines of the pressure around the wing 
at several span stations. 



Figure 11. Magnified deformations and pressure 
contour lines at several stations of the wing. Two 
frames during the oscillation are shown. 

The oscillations of the wing tip are shown in figure 12 
for two different values of the dynamic pressure normal- 
ized with the experimentally observed flutter value. The 
figure demonstrates that the predicted oscillations are 
damped at 90% of the experimental flutter pressure and 
grow at 1 10% of that value. 



Figure 12. Time history of the tip displacement for 
the wing at Mach number 0.9 using two values of 
the dynamic pressure. 


8,3, Impeller Flow 

The underlying function of turbomachinery equipment, 
such as pumps, compressors, turbines or fans, is to 
smoothly impart, or extract, energy from rotating blades 
to increase, or decrease, the velocity and pressure of a 
fluid stream. Proper aerodynamic and structural blade 
design is therefore critical for achieving optimal perfor- 
mance. 



Figure 13. Geometry of the centrifugal impeller 
study. 

The utility of multi-disciplinary analysis and optimiza- 
tion tools for design applications in aircraft propulsion 
is currently being investigated by the NASA Lewis 
Research Center. One of the turbomachinery applica- 
tions involves the simulation of fluid flow through the 
rotor ot an Allied Signal centrifugal compressor. A vis- 
cous, steady, compressible flow finite element analysis 
was performed, in which turbulence was modeled using 
the Spalart-Allmaras model. 

The gometry of the centrifugal impeller analyzed in 
this study is shown in figure 13. It consists of 17 main 
and 17 splitter blades located with equal spacing along 
the circumferential direction. The axial length of the 
impeller is 3.64 in., whereas its inlet and outer diameters 
are 2.20 in. and 8.4 in., respectively. The clearance 
between the main blades and the shroud surface is 
0.0033 in. at the leading edge of the blades and 0.0120 
in. at the trailing edge of the blades. 

The relative frame of reference employed in the current 
analysis is taken with respect to an observer sitting on 
the hub and rotating fixed with the wheel around the 
axial axis of the compressor. In this relative frame of 
referent e, the hub, main and splitter blades are station- 
ary anc the shroud is instead rotating in the reverse 
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direction. This approach greatly simplifies the numerical 
analysis, because it allows the modeling of the fluid flow 
passage on a grid that does not change with time, since 
the shroud has an ax i symmetric geometry with respect 
to the axis of rotation. In addition, due to the cyclic sym- 
metry of the flow, the modeling of only one rotor (main) 
blade passage is sufficient, which greatly reduces CPU 
time and memory requirements. For clarification, the 
notation “cyclic symmetry” in a Cartesian coordinate 
system corresponds to “periodicity” in a cylindrical 
coordinate system. 


Figure 14. Pressure contours on the flow domain of 
the centrifugal impeller. 

The pressure distribution on several k-surfaces (cross 
flow surfaces) is shown in figure 14. The pressure rises 
about 3 atm. in a uniform fashion from the inlet to the 
exit. As can be seen in the close-up, the pressure is 
indeed lower on the suction side of the main blade and 
splitter blade and increases at the pressure side. Further 
examination revealed small loading on the blade. This is 
expected, since an off-design condition was considered 
here. No unloading, however, was found other than at 
the trailing edge region of the blade. 

Contours of the relative Mach number on several blade- 
to-blade surfaces are shown in figure 15. For the case 
considered here, the incoming flow at the inlet is at sub- 
sonic conditions. The high relative Mach number values 
seen at the inlet close to the shroud are due to the contri- 
bution from the high relative tangential velocity compo- 
nent. The flow appears to become locally transonic at 
the leading edge of the main blades, which is most likely 


due to the modeling of the leading edge as a sharp cor- 
ner rather than as a smooth elliptical surface. Transonic 
regions are also found near the hub and between the 
main and splitter blades due to the presence of oblique 
shock waves. Finally, a recirculation region was found 
that originated upstream of the leading edge of the main 
blade close to the shroud. 



Figure 15. Mach number contours at different 
stations within the flow domain of the centrifugal 
impeller. 

8.4. Hypersonic Cowl Lip Cooling 

The analysis of high aerothermodynamic loading on the 
leading edges of future hypersonic flight vehicles leads 
to a clear thermal FSI problem. For example, to prevent 
material failure on the engine cowl lip caused by aero- 
dynamic heating (including the shock-on-shock interac- 
tion) that occurs at hypersonic flight range, active 
internal cooling of the leading edges is deemed neces- 
sary. As a result, this problem class encompasses exter- 
nal hypersonic (compressible) flow with high heating 
rates near stagnation points, internal coolant (incom- 
pressible) flow through the leading edge with heat con- 
duction and convection, and thermal conduction and 
thermal stresses in the leading edge structure. 

Melis et.al. 17 have conducted experiments in the NASA 
Lewis Hot Gas Facility to evaluate alternate design con- 
cepts of internal coolant passage for the cowl lip. A sin- 
gle-physics analysis to compute the material thermal 
stress, whereby the external aerothermodynamic loading 
and internal coolant heat transfer were approximated by 
empirical correlations, has also been performed. Results 
were found to be encouraging. Nevertheless, a design 
tool which incorporates direct multiphysics analysis is 
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still warranted to accurately characterize the leading 
edge models. 



Figure 16. Actively cooled leading edge geometry 
model. 

Work is currently underway to employ Spectrum to per- 
form multiphysics analysis on the impingement-cooled 
cowl lip model made of nickel. Figure 16 depicts the 
actively cooled leading edge geometry. The internal 
coolant (gaseous hydrogen) impinges on the leading 
edge from the center passage, carries away the heat flux, 
and returns to the outflow manifold along the lip sur- 
face. The coolant is running at high speed but still falls 
in the incompressible turbulent flow range. 




Figure 17. Solid and internal coolant grids for the 
impingement cooled cowl lip model. 

The mesh discretizations for the internal turbulent cool- 
ant fluid and the solid component in this fluid-solid 
interaction problem is shown in Figure 17. Note that the 
mesh for the internal coolant is highly stretched near the 
channel walls to account for the sharp gradient in the 
turbulent boundary layer 

Initial results for the steady-state coolant velocity vector 
fields and the temperature distribution of the solid are 
shown in figures 18 and 19. To validate the current anal- 
ysis, external compressible flow interaction was not 
considered first. Instead, the external heat loading was 


simulated by imposing spatial heat flux on the solid sur- 
face. It is seen that the flow field has been well resolved, 
especially the stagnation region on the tip and the comer 
recirculation region where the flow changes direction 
are clearly identified. Although turbulence results are 
not shown, the wall turbulent boundary layer has also 
been resolved adequately. 



Figure 18. Velocity vectors for the internal gaseous 
hydrogen flow of the leading edge cooling. 



Figure 19. Temperature contours of the solid 
(nickel) in the leading edge cooling. 


In Figure 19, the temperature distribution is found to 
exhibit a strong gradient in the nose region where the 
external heat loading is at the maximum. The maximum 
tip temperature is about 1300 degrees R and levels off to 
600 decrees downstream. The inlet coolant temperature 
is set at 500 degrees and it reaches a maximum tempera- 
ture of {.bout 750 degrees in the nose region. 

8,5. Quenching Process Simulation 

A quenching process is relatively rapid cooling of a 
workpiece by immersion into a fluid after some high- 
temperature metal-forming operation. It is an integral 
part of i manufacturing process as it helps determine the 
materia properties of the finished component created 
from the.: workpiece. Quenching is a multi-physics prob- 
lem: heat transfer from the solid to the fluid induces 
flow which alters the cooling behavior. Moreover, the 
overall objective of a complete simulation would be to 
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include ihe effects of solid deformation and nonlinear 
material response such as plasticity. 



Figure 20. Cross-sectional dimensions of a steel 
forging (axisymmetric). 

For this initial investigation, a thermal-only quench 
problem for a simple disk was specified. The physical 
dimensions of the forging are indicated in figure 20. The 
overall geometry of the quench tank and the location of 
the immersed forging are illustrated in figure 21. The 
idealized simulation begins with the disk already in the 
tank at a temperature of 2000 °F surrounded by fluid at 
107 °F. The subsequent five minutes of cooling are 
assumed to ignore heat transfer into the air from the oil 
surface. Effects associated with boiling of the quenchant 
were also neglected. 


I 

I 



9 in 


Figure 21. Schematic cross-section of the quench 
tank. 

The quenchant is modeled as an incompressible fluid. 
The effect of thermal buoyancy is captured by the 
Boussinesq approximation, which maintains a constant 
density assumption in the equations of motion while 
adding an effective body force proportional to the 
expansivity of the fluid and the local temperature varia- 
tion from the reference (initial) value. The surfaces of 


the fluid domain corresponding to the tank walls arc 
designated with no-slip velocity boundary conditions 
and held at a constant temperature equal to the initial 
fluid temperature of 107 °F. Free surface effects were 
not included. Instead, the top surface was considered 
fixed with the vertical velocity set to zero. This surface 
was considered adiabatic, so no effects of conduction/ 
convection into the surrounding air are modeled. Tem- 
perature dependant properties are used to model the oil 
and steel materials. Figure 22 presents a closc-up of the 
mesh in the vicinity of the forging. 



Figure 22. Close up of forging mesh. 

In this simulation, automatic time incrementation with a 
maximum step size of 1 .25 seconds was used to model a 
300 second forging interval. Figures 23-25 display the 
system response at time 20.2 seconds. Figure 23 depicts 
contours of temperature (°F) in the forging and quen- 
chant. Figure 24 shows contours of the magnitude (ft / 
sec) of the vertical velocity in the quenchant. 

The early time response shows a thermal plume devel- 
oping from the rim of the forging due to heat transfer 
along the vertical surfaces. This ring-like structure then 
collapses into a central plume, which is likely to be at 
least in part an artifact of the axisymmetric assumption 
of the simulation. Multiple plumes rise from the top hor- 
izontal surface and are also subsumed into this central 
plume. The central plume rises to the surface setting up 
a recirculation cell in the upper half of the tank. The 
healed fluid at the bottom horizontal surface of the forg- 
ing forms a thinner boundary layer of low velocity as 
the induced buoyancy “traps” the liquid against that sur- 
face. It is at this location that neglecting phase change 
(boiling) of the quenchant is likely most detrimental. 
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Figure 23, Temperature contours in the forging and 
quenchant at time 20.2 seconds. 




after 20.2 seconds of quench. The outboard edges of the 
forging cool most rapidly, as seems intuitive. Perhaps 
most interesting is some of the local structure exhibited 
during early lime in the temperature contours near the 
top surface. Here local variations in the quenchant 
velocities lead to differing cooling rates and local cool 
spots. Such a mechanism could lead to local variations 
in material properties which might be a performance 
issue for the final machined part. 
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Figure 24. Vertical fluid velocity contours in the 
quenchant at time 20.2 seconds. 

The overall objective of a complete quenching analysis 
is to determine the thermal history of the forging and 
hence to subsequently predict the mechanical properties 
“locked in” the material after the cooling. Figure 25 
shows in isolation the thermal contours in the workpiece 


Figure 25. Temperature contours in the forging at 
lime 20.2 seconds. 

Conclusions 

A multiphysics simulation approach based on the finite 
element method has been described. This work 
addresses compressible and incompressible fluid flow, 
structural, and thermal modeling as well as the interac- 
tion between these disciplines. The approach is based on 
a single computational framework for the modeling of 
multiple interacting physical phenomena. The aug- 
mented-Lagrangian method is used to enforce interac- 
tion constraints among all field variables in a fully- 
coupled manner. Consistent finite element treatments of 
uniforn region balance laws were described within the 
multiply sics framework. The arbitrary-Lagrangian- 
Eulerian method is utilized to account for deformable 
fluid domains. 

The efficacy of this method in simulating coupled fluid- 
solid-thermal interaction was demonstrated with 
aeroelastic, thermal management and flow-induced 
vibration problems. The excellent scalability of this 
approa:h was illustrated. 
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